Glacier mass balance

We simplify the glacier mass balance to the following equation

\[ \Delta S = P-M = A_{imbal} \]

where the change of water storage (\(\Delta S\)) is equal to the precipitation (\(P\), assumed solid) minus the glacier melt (\(M\)). Typically melt excesses precipitation and we have negative imbalance ablation, indicating glacier storage loss. The glacier mass balance is calculated in annual time steps. We thereby refer to the hydrological year starting on October 1st of the previous year to take a full accumulation and ablation season into account.

Temperature index model

Hock (2003) describes several variants of the temperature index model for simulating glacier melt. The package implements the temperature index melt model described in (Hock 2003) in the function .

\[ M = \biggl\{ \begin{array}{l, l} 0, & T < T_{threshold} \\ f_{M} \cdot \left( T - T_{threshold} \right), & T >= T_{threshold} \end{array} \]

where \(M\) is the glacier melt in \(mm/d\), \(T\) is the daily average temperature in \(^{\circ} C\). The two parameters \(f_{M}\) and \(T_{\text{threshold}}\) refer to the melt factor and the threshold temperature above which glacier melt occurs and need to be calibrated. They have the units \(\frac{mm}{^{\circ} C \cdot d}\) and \(^{\circ} C\) respectively. Glacier melt is calculated in daily time steps.

Glacier volume development

As glaciers melt, their volume changes. This has to be taken into account for the long-term simulation of glacier discharge. To determine the initial glacier volume we use the area of the geometry of the RGI v6.0 data set and multiply it with the average thickness of the glacier. Glaciers larger than 1km2 we sub-divide into elevation bands of 100m altitude to account for elevation dependent temperature forcing. The large glaciers melt from the lowest elevation band to the highest elevation band whereby the glacier melt is subtracted from the glacier volume of lowest elevation band which is still glacierized. The small glaciers are not spatially discretized and thus they are melted homogeneously.

For annual time step \(t\), the evolution of the glacier volume is calculated as follows:

\[ A(t) = \text{glacierArea_RGIF}\bigl(V(t)\bigr) \]

\[ Q_{glacier}(t) = M(t) \cdot A(t) \]

\[ V(t+1) = V(t) + \Delta S = V(t) + \text{glacierImbalAbl}\bigl(M(t)\bigr) \]

Note that the imbalance ablation (or \(\Delta S\)) is negative for glacier mass loss and positive for glacier accumulation.

\(Q_{\text{glacier}}\) can be calibrated against glacier discharge derived from the Miles & Hugonnet data sets. The automated calibration is not yet included in the package.

The function is an empirical scaling function analogue the inverse of the scaling function derived by Erasov (1968) but based on the modern RGI v6.0 glacier geometries and the glacier thickness data set by Farinotti et al. (2019). The package implements volume-area and area-volume scaling functions based on both, Erashov and RGI-Farinotti data, allowing the estimation of glacier areas based on glacier volumes ( and ) and estimations of glacier volumes based on glacier areas ( and ).

The function is an empirical scaling function relating glacier imbalance ablation, i.e. the glacier storage loss glacier melt. It is derived from the glacier discharge data set by Miles et al. (2021) and the glacier thinning data set by Hugonnet et al. (2021).

Propagation of uncertainty

The initial glacier volume of each glacier is attributed an uncertainty of p/m 26%. This number is based on the average uncertainty of the glacier volume per RGI region reported in Farinotti et al. (2019) (Table 1). The error is the same for RGI regions 13 and 14. The uncertainty of the area of the RGI v6.0 glacier outlines is not known. We therefore assume that the error of the glacier volume stems to 50% from the estimation of the glacier thickness and to 50% from the glacier area. We further assume that the errors of the glacier area and glacier thickness are un-correlated and can thus estimate the uncertainties of the glacier area data and the glacier thickness data to be p/m 13% each.

For the non-linear relationships and , the standard deviation of the residuals of the fit was computed. The estimated error of the fit is assumed to be equal plus/minus twice the standard deviation of the residuals and yields 31% and 53% respectively. The residuals are not normally distributed and their actual distribution in unknown. Nevertheless, we apply the principles of linear error propagation for the sake of simplicity. We assume uncorrelated errors and therefore add the error of the function input to the error of the fit to obtain the error of the function output.

\[ \varepsilon_{V} = \varepsilon_{A} + \varepsilon_{\text{glacierVolume_RGIF}} = 0.26 + 0.31 = 0.57 \]

\[ \varepsilon_{A} = \varepsilon_{V} + \varepsilon_{\text{glacierArea_RGIF}} = 0.26 + 0.53 = 0.79 \]

We do not have an estimate for the error of the temperature index model. We assume a conservative relative error of 2, indicating that the estimated glacier melt is within a range of plus/minus 2 times it’s value.

The error of the imbalance ablation amounts to 73%. Assuming independent errors from the temperature index model for glacier melt and the scaling function between imbalance ablation and glacier melt, the error of the imbalance ablation amounts to approximately:

\[ \varepsilon_{Q_{\text{imb,melt}}} = \varepsilon_{Q_{\text{M}}} + \varepsilon_{\text{glacierImbalAbl}} = 2 + 0.73 = 2.73 \]

The errors stated above relate to initial estimates of the glacier areas, volumes, total ablation and imbalance ablation. The uncertainties are not propagated through each time step of the model as the errors depend on the errors of the previous time step, i.e. they are not uncorrelated over time. In any case, the uncertainty margins for any glacier melt modelling done with the presented method are large.

# Relative error estimates for the initial state in %
error_stats = tibble::tibble(
  sV = 0.26,  # Glacier volume
  sA = 0.13,  # Glacier area
  sTh = 0.13,  # Glacier thickness
  sMelt = 2,  # Glacier melt
  sVolumeRGIF = 0.31,  # Calculating glacier volume from glacier area
  sAreaRGIF = 0.53,  # Calculating glacier area from glacier volume
  sImbal = 0.73,  # Calculating glacier imbalance ablation from glacier melt
  sQglgrowth = 1.6,  # Calculating total ablation from glacier thinning  
  sQglmelt = 1.6)  

Demonstration

We demonstrate the above described method with the data from the Atbashy basin, a tributary of the Naryn river. As for the previous chapter you can download all required data from this link.

library(tmap)
library(timetk)
library(sf)
library(raster)
library(exactextractr)
library(tidyverse)
library(lubridate)

devtools::install_github("hydrosolutions/riversCentralAsia")
library(riversCentralAsia)

# Stochastic optimisation algorithm for model calibration
library(GA)

# Path to the data directory downloaded from the download link provided above. 
# Here the data is extracted to a folder called atbashy_glacier_demo_data
data_path <- "../../atbashy_glacier_demo_data/"

Spatial discretization

Glacier melt is calculated on elevation bands. The following code chunks show how to derive elevation bands on glaciers.

Glaciers smaller than 1 km2 are not sub-divided to limit the use of computational resources. The area threshold of 1 km2 is arbitrary but similar area thresholds have been used in the literature (e.g. Miles et al. (2021)).

library(sf)
library(raster)
library(exactextractr)
library(tidyverse)
devtools::install_github("hydrosolutions/riversCentralAsia")
library(riversCentralAsia)

data_path <- "../../atbashy_glacier_demo_data/"

# Define output file name
outfile <- paste0(data_path, "GIS/rgi_glaciers_atbaschy_el_bands.shp")

# Read in HRUs for the HBV models
hru <- st_read(paste0(data_path, "GIS/16076_HRU.shp"))

# Read in glacier geometries. Here we have already intersected the glacier 
# outlines from RGI v6.0 with the outlines of the sub-basins of the Atbashy 
# model in QGIS, so we know which glacier lies in which sub-basin. For the 
# generation of the elevation bands, you can also use a simple RGI layer. 
rgi <- st_read(paste0(data_path, "GIS/16076_Glaciers_per_subbasin.shp")) |> 
  dplyr::select(RGIId, Area, Zmed) # Area in km2, Zmed in masl (WGS84)

rgi_large_glaciers <- rgi |> filter(Area >= 1) 
rgi_small_glaciers <- rgi |> filter(Area < 1) 

dem <- raster(paste0(data_path, "GIS/16076_DEM.tif")) |>
  mask(rgi_large_glaciers) |> 
  setMinMax()

rgi_large_glaciers <- rgi_large_glaciers |> st_transform(crs = 4326)
rgi_small_glaciers <- rgi_small_glaciers |> st_transform(crs = 4326)


## Elevation Band Zonation Parameters
band_interval <- 100 # in meters
holeSize_km2 <- 0.2    # cleaning holes smaller than that size
smoothFact <- .5     # level of band smoothing
demAggFact <- 1      # dem aggregation factor (carefully fine-tune this)

## Generate elevation bands for the large glaciers
hru_shp <- gen_basinElevationBands(
  dem_PathN = dem,
  demAggFact = demAggFact,
  band_interval = band_interval,
  holeSize_km2 = holeSize_km2,
  smoothFact = smoothFact)
hru_shp_latlon <- hru_shp |> 
  st_transform(crs = 4326) |> 
  st_make_valid() # now latlon

## Intersect the elevation bands with the rgi layer
rgi_large_glaciers_el_bands <- st_intersection(rgi_large_glaciers, hru_shp_latlon)

## Concatenate with the small glaciers
rgi_elbands <- rgi_large_glaciers_el_bands |> 
  add_row(rgi_small_glaciers |> 
            mutate(layer = 1))

# Now we have a GIS layer with small glacier represented by 1 elevation band and 
# large glaciers represented by several elevation bands
st_write(hru_out, outfile, append = FALSE)

We can do some post-processing of the GIS layer to add an estimate of the glacier volume and the areas of the elevation bands.

# Read DEM 
dem <- raster(paste0(data_path, "GIS/16076_DEM.tif"))

# Read in HRUs for the HBV models
hru <- st_read(paste0(data_path, "GIS/16076_HRU.shp"))
## Reading layer `16076_HRU' from data source 
##   `C:\Users\marti\hydrosolutions Dropbox\Bea martibeatrice@gmail.com\marti\github\atbashy_glacier_demo_data\GIS\16076_HRU.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 40 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1105700 ymin: 4575416 xmax: 1197076 ymax: 4616616
## Projected CRS: WGS 84 / UTM zone 42N
# Read in glacier geometries. Here we have already intersected the glacier 
# outlines from RGI v6.0 with the outlines of the sub-basins of the Atbashy 
# model in QGIS, so we know which glacier lies in which sub-basin.  
rgi <- st_read(paste0(data_path, "GIS/16076_Glaciers_per_subbasin.shp")) |> 
  # Area in km2, Zmed in masl (WGS84)
  dplyr::select(RGIId, Area, Zmed, Zmax, Slope, Aspect) |> 
  distinct(RGIId, .keep_all = TRUE)
## Reading layer `16076_Glaciers_per_subbasin' from data source 
##   `C:\Users\marti\hydrosolutions Dropbox\Bea martibeatrice@gmail.com\marti\github\atbashy_glacier_demo_data\GIS\16076_Glaciers_per_subbasin.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 135 features and 24 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 76.2421 ymin: 41.11981 xmax: 77.3168 ymax: 41.41837
## Geodetic CRS:  WGS 84
# Read in the glacier elevation bands layer
rgi_elbands <- st_read(paste0(data_path, 
                              "GIS/rgi_glaciers_atbaschy_el_bands.shp"), 
               quiet = TRUE) |> 
  st_transform(crs = crs(dem)) 

# Load a pre-processed raster file with glacier thickness in the Atbashy basin
# (see vignette [glaciers 01]{glaciers-01-intro.html} for details on the pre
# -processing)
glacier_thickness <- raster(
  paste0(data_path, 
  "GLACIERS/Farinotti/pre-processed_glacier_thickness.tif"))

# Only keep the variables we need for this analysis
rgi_elbands <- rgi_elbands |> 
  dplyr::select(RGIId, Area, elvtn_b) |> 
  rename(Area_tot_glacier_km2 = Area) |> 
  mutate(ID = paste0(RGIId, "_", elvtn_b)) 

# Get mean elevation of each glacier/elevation band from DEM
rgi_elbands$z_masl <- exact_extract(dem, rgi_elbands, "mean", progress = FALSE)

# Update the glacier area within the basin boundaries
rgi_elbands$A_km2 <- as.numeric(st_area(rgi_elbands))*10^(-6)
glaciers <- unique(rgi_elbands$RGIId)

# Calculate the average glacier thickness for each elevation band. 
rgi_elbands$thickness_m = exact_extract(glacier_thickness, 
                                        rgi_elbands, "mean", progress = FALSE)
rgi_elbands$thickness_m[is.na(rgi_elbands$thickness_m)] <- 0

tmap_mode("view")
tm_shape(dem, name = "DEM") + 
  tm_raster(n = 6, 
            palette = terrain.colors(6),
            alpha = 0.8,
            legend.show = TRUE, 
            title = "Elevation (masl)") + 
  tm_shape(rgi_elbands, 
           name = "RGI v6.0 in elevation bands") + 
  tm_polygons(col = "z_masl", lwd = 0.4, 
              title = "Mean glacier\nelevation [masl]", 
              palette = "Blues") + 
  tm_scale_bar(position = c("left", "bottom"))

Observed glacier melt

Observed average glacier thinning rates by Hugonnet et al. (2021) are shown in the figure below. 30 out of 135 glaciers in the Atbashy basin show increase of glacier thickness between 2000 and 2020.

# Load the glacier thinning data set and filter it to the glaciers in the 
# catchment of interest.  
hugonnet <- read_csv(paste0(
  data_path, "/GLACIERS/Hugonnet/dh_13_rgi60_pergla_rates.csv"), 
  show_col_types = FALSE)
hugonnet <- hugonnet |> 
  dplyr::filter(rgiid %in% rgi$RGIId) |> 
  tidyr::separate(period, c("start", "end"), sep = "_") |> 
  mutate(start = as_date(start, format = "%Y-%m-%d"), 
         end = as_date(end, format = "%Y-%m-%d"), 
         period = round(as.numeric(end - start, units = "days")/366))
glaciers_hugonnet <- rgi |> 
  left_join(hugonnet |> dplyr::select(rgiid, area, start, end, dhdt, err_dhdt, 
                                      dvoldt, err_dvoldt, dmdt, err_dmdt, 
                                      dmdtda, err_dmdtda, period),  
            by = c("RGIId" = "rgiid")) 

# Add glacier discharge by Miles et al.
miles <- read_delim(paste0(data_path,
                           "GLACIERS/Miles/Miles2021_Glaciers_summarytable_20210721.csv"),
                    show_col_types = FALSE, delim = ",") |> 
  dplyr::filter(RGIID %in% rgi$RGIId & VALID == 1)

glaciers_hugonnet <- glaciers_hugonnet |> 
  left_join(miles |> dplyr::select(RGIID, totAbl, totAblsig, imbalAbl, 
                                   imbalAblsig), 
            by = c("RGIId" = "RGIID")) |> 
  mutate(Qgl_m3a = ifelse(is.na(totAbl), glacierDischarge_HM(dhdt), -totAbl))

observed_melt <- glaciers_hugonnet |> 
  st_drop_geometry() |>
  dplyr::filter(period == 1) |> 
  dplyr::select(RGIId, start, dhdt, Qgl_m3a, Area) |> 
  mutate(Area_m2 = Area * 10^6) |> 
  dplyr::select(-Area) |> 
  distinct()

totQgl <- glaciers_hugonnet |> 
  st_drop_geometry() |> 
  dplyr::filter(period == 20) |> 
  mutate(Qgl_m3a = ifelse(Qgl_m3a < 0, 0, Qgl_m3a), 
         dVweqda_m3a = ifelse(dmdtda<0, -dmdtda, 0) * area) |> 
  summarise(Qgl_m3a = sum(Qgl_m3a), 
            dVweqda_m3a = sum(dVweqda_m3a)) |> 
  transmute(Qgl_mcma = Qgl_m3a * 10^(-6), 
            Qgl_m3s = Qgl_m3a / (60*60*24*365), 
            dVweqda_mcma = dVweqda_m3a * 10^(-6), 
            dVweqda_m3s = dVweqda_m3a / (60*60*24*365))

tmap_mode("view")
tm_shape(dem, name = "DEM") + 
  tm_raster(n = 6, 
            palette = terrain.colors(6),
            alpha = 0.8,
            legend.show = TRUE, 
            title = "Elevation (masl)") + 
  tm_shape(glaciers_hugonnet |> dplyr::filter(period == 20) |> 
             mutate(Qgl_mcma = Qgl_m3a * 10^(-6)), 
           name = "Glacier discharge\n[10^6 m3/a]") + 
  tm_polygons(col = "Qgl_mcma", lwd = 0.4, 
              title = "Qgl (10^6 m3/a)") + 
  tm_scale_bar(position = c("right", "bottom"))

The average annual contribution of glacier melt to river runoff is 1.9198103 m3/s. The annual norm discharge at Atbashy gauging station is 16.6 m3/s which yields an average annual glacier contribution of about 11.5651221 % which is reasonable given the high median altitude of the basin of 3470 masl.

Forcing

There will be a separate Vignette to demonstrate how to prepare the forcing. For now, please refer to the methodology described in the CAHAM book. We provide you with the pre-processed forcing data and give you just a brief overview over the data source. We have sub-divided the glaciers into elevation bands to be able to reproduce the elevation-dependent temperature forcing on the glacier.

Historical temperatures

As meteorological data for high elevations in Central Asia is very scarce we use the CHELSA v2.1 data set (Karger et al. 2017). This is a global data set of forcing data for hydrolgical models based on ERA5 but corrected for precipitation biases in mountainous regions and especially suited for high elevations. The daily CHELSA forcing has been cut to the Central Region by the originator of the data set, D. Karger, WSL and extracted to the hydrological response units of the glaciers in the Atbashy basin by T. Siegfried.

hist_obs <- readRDS(file = paste0(data_path, "CLIMATE/hist_obs_glacier_tas.rds"))

# Plot the temperature time series for a given glacier/elevation band
hist_obs[, c(1, 2, 5, 10)] |> 
  pivot_longer(-date, names_to = "layer", values_to = "Temp_degC") |>
  mutate(layer = as.numeric(gsub(pattern = "_", "", 
                                 str_extract(layer, "_\\d$")))) |> 
  plot_time_series(
    .date_var = date, 
    .value = Temp_degC, 
    .color_var = layer, 
    .color_lab = "El. band", 
    .plotly_slider = TRUE, 
    .smooth = FALSE, 
    .line_size = 0.5
  )

Daily temperature time series of the lowest, a middle-elevation and the highest elevation band of glacier RGI60-13.08930 in the Atbashy basin. Data source: CHELSA v2.1.

Please note that, as expected, the daily temperatures are higher in lower elevation bands (lower elevation band numbers) than in higher elevation bands.

Future temperatures

Future temperature development per glacier or elevation band is extracted from the 3 CMIP6 GCM models with highest priorities for the region (Krasting et al. 2018; Yukimoto et al. 2019; Tang et al. 2019) downloaded from COPERNICUS. We take 4 socioeconomic scenarios (SSP) into account, covering 4 different emission scenarios (RCP): SSP 1 for RCP 2.6, short ssp126; SSP 2 for RCP 4.5, short ssp245; SSP 3 for RCP 7.0, short ssp370; SSP 5 for RCP 8.5, short ssp585. The temperatures of the climate models are bias corrected using the CHELSA data and a quantile mapping method (Gudmundsson et al. 2012; Gudmundsson 2016). The figure below shows the development of the temperature in the lowest elevation band of glacier RGI60-13.08930 as an example.

fut_sim <- readRDS(file = paste0(data_path, "CLIMATE/fut_sim_glacier_tas_qmapped.RDS"))

# Plot the temperature time series for a given glacier/elevation band
glacier <- "RGI60-13.08930_1"

# Extract the temperature for the selected glaciers for all GCMs and SSPs
scenarios <- names(fut_sim)
fut_temp <- NULL
for (scenario in scenarios) {
  fut_temp <- rbind(fut_temp, 
                    fut_sim[[scenario]] |> 
                      dplyr::select(Date, all_of(glacier)) |> 
                      mutate(Scenario = scenario))
}
fut_temp <- fut_temp |> 
  mutate(Hyear = hyear(Date)) |> 
  group_by(Hyear, Scenario) |> 
  summarise(Date = first(Date), 
            Temp = mean(get(glacier))) |> 
  separate(Scenario, into = c("GCM", "SSP"), sep = "_") |> 
  ungroup() |>
  dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear), 
                GCM != "IPSL-CM6A-LR")

# Plot annual data for better readability
ggplot() + 
  geom_line(data = hist_obs |> 
              mutate(Hyear = hyear(date)) |> 
              group_by(Hyear) |> 
              summarise(date = first(date), 
                        Temp = mean(get(glacier))) |> 
              ungroup() |> 
              dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)), 
            aes(date, Temp)) + 
  geom_line(data = fut_temp, aes(Date, Temp, colour = GCM, 
            linetype = SSP)) + 
  scale_colour_viridis_d() + 
  labs(x = "Date", y = "T [deg C]", title = glacier) +
  theme_bw()
Annual historical (CHELSA) and future (CMIP6, Copernicus) temperature for the lowest elevation band of glacier RGI60-13.08930.

Annual historical (CHELSA) and future (CMIP6, Copernicus) temperature for the lowest elevation band of glacier RGI60-13.08930.

Calculate daily glacier melt

The code snippet below illustrates how to calculate the daily glacier melt in a river catchment with default parameters for the temperature index model. The figure shows daily melt rates for the different elevation bands of one of the larger glaciers (RGI60-13.08930). The highest melt occurs at low elevation bands. Glacier melt only occurs during the summer months when the average daily temperature is above the threshold temperature for glacier melt.

# Default parameters of the temperature index model for glacier melt. 
MF_small = 2  # Melt coefficient for small glaciers (smaller than 1km2)
MF_large = 0.1  # Melt coefficient for large glaciers (larger or equal than 1km2)
threshold_temperature = 0  # Default temperature threshold for glacier melt 

Area <- rgi_elbands |> 
  st_drop_geometry() |> 
  dplyr::select(ID, A_km2) |> 
  pivot_wider(names_from = ID, values_from = A_km2)
# Assign different melt factors to large and small glaciers. 
MF <- Area |> 
  mutate(across(everything(), ~MF_large), 
         across(ends_with("_1"), ~MF_small))

melt <- glacierMelt_TI(temperature = hist_obs |> dplyr::select(-date),
                       MF = MF,
                       threshold_temperature = threshold_temperature)
melt <- as_tibble(melt) |> 
  mutate(date = hist_obs$date) |> 
  relocate(date, .before =  where(is.numeric))

ggplot(melt[, 1:9] |> 
         pivot_longer(-date, names_to = "ID", values_to = "Melt") |> 
         left_join(rgi_elbands |> 
                     st_drop_geometry() |> 
                     dplyr::select(ID, z_masl), by = "ID") |> 
         mutate("Mean\nelevation\n[masl]" = as.factor(round(z_masl)))) + 
  geom_line(aes(date, Melt, colour = `Mean\nelevation\n[masl]`)) + 
  scale_colour_viridis_d() + 
  labs(x = "Date", y = "Melt [mm/d]") + 
  theme_bw()
Daily glacier melt per elevation band.

Daily glacier melt per elevation band.

Compare to measured glacier melt

Aggregate the daily glacier melt to annual data and compare it to the observed glacier melt.

melt_a_eb <- melt |> 
  pivot_longer(-date, names_to = "ID", values_to = "Melt") |> 
  mutate(Hyear = hyear(date)) |> 
  group_by(Hyear, ID) |> 
  summarise(date = first(date), 
            Melt = sum(Melt), 
            .lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0, 
                              Melt*(1-error_stats$sMelt)), 
            .ub_Melt = Melt*(1+error_stats$sMelt))|> 
  ungroup() 
melt_a <- melt_a_eb |> 
  separate(ID, into = c("RGIId", "elB"), sep = "_") |> 
  group_by(Hyear, RGIId) |> 
  summarise(date = as_date(first(date)), 
            Melt = sum(Melt), 
            .lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0, 
                              Melt*(1-error_stats$sMelt)), 
            .ub_Melt = Melt*(1+error_stats$sMelt)) |> 
  ungroup() |> 
  dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)) |> 
  left_join(rgi |> dplyr::select(RGIId, Area) |> st_drop_geometry(), 
            by = "RGIId") |> 
  # Assume constant glacier area over the comparison period. 
  mutate(Qgl_m3a = Melt*10^(-3) * Area*10^6, 
         .lb_Qgl_m3a = Qgl_m3a * (1-sqrt(error_stats$sMelt^2+error_stats$sA^2)),   
         .lb_Qgl_m3a = ifelse(.lb_Qgl_m3a < 0, 0, .lb_Qgl_m3a), 
         .ub_Qgl_m3a = Qgl_m3a * (1+sqrt(error_stats$sMelt^2+error_stats$sA^2)))

melt_obs_a_all <- observed_melt |> 
  mutate(.lb_Qgl_m3a = ifelse(
           dhdt > 0, 
           ifelse(Qgl_m3a*(1-error_stats$sQglgrowth)<0, 0, 
                  Qgl_m3a*(1-error_stats$sQglgrowth)), 
           ifelse(Qgl_m3a*(1-error_stats$sQglmelt)<0, 0, 
                  Qgl_m3a*(1-error_stats$sQglmelt))),
         .ub_Qgl_m3a = ifelse(dhdt > 0, 
                              Qgl_m3a*(1+error_stats$sQglgrowth), 
                              Qgl_m3a*(1+error_stats$sQglmelt)))
melt_obs_a <- melt_obs_a_all |> 
  dplyr::filter(RGIId %in% glaciers[c(6, 7)])
  
ggplot() + 
  geom_ribbon(data = melt_a |> 
                dplyr::filter(RGIId %in% glaciers[c(6,7)]), 
              aes(date, Qgl_m3a, ymin = .lb_Qgl_m3a, ymax = .ub_Qgl_m3a, 
                  fill = RGIId), colour = NA, alpha = 0.1) + 
  geom_ribbon(data = melt_obs_a, 
              aes(start, Qgl_m3a, 
                  ymin = .lb_Qgl_m3a, 
                  ymax = .ub_Qgl_m3a, 
                  colour = RGIId, linetype = "obs", fill = RGIId), 
              size = 0.2, alpha = 0.1) + 
  geom_line(data = melt_a |> 
              dplyr::filter(RGIId %in% glaciers[c(6,7)]), 
            aes(date, Qgl_m3a, colour = RGIId, linetype = "sim")) + 
  geom_line(data = melt_obs_a, aes(start, Qgl_m3a, colour = RGIId, 
                                   linetype = "obs")) + 
  labs(x = "Date", y = "Melt [m weq/a]") + 
  scale_linetype_manual(name = "Source", 
                        values = c("sim" = 1, "obs" = 2)) + 
  scale_colour_viridis_d() + 
  scale_fill_viridis_d() + 
  ggtitle(paste0("MF: ", MF, "Tth: ", threshold_temperature)) + 
  theme_bw()
Daily glacier melt per elevation band.

Daily glacier melt per elevation band.

The temperature index model can be calibrated with the specific glacier volume change provided by Hugonnet et al. (2021) (see \(\ref{fig:glaciermeltcomparison}\)).

Calibration of the temperature index model

The 2 parameters of the temperature index model can be calibrated using observed glacier melt rates. As mentioned above, the observed glacier melt rates are simulated total ablation by Miles et al. (2021) where available or estimated total ablation using an empirical relationship between total glacier ablation by Miles et al. (2021) and glacier thinning rates by Hugonnet et al. (2021).

The calibration is done using an R implementation of a genetic algorithm Scrucca (2017). The optimization is slow but it is necessary to use a global optimizer (nlminb of the stats package for example will not converge for all glaciers).

# Calibrating many glaciers takes a lot of time. 
# Reduce the amount of time required and only calibrate the glaciers that are 
# larger or equal than 0.2 km2. That represents 83% of the glacier volume in the 
# Atbashy basin but reduces the number of glaciers to be calibrated by 65%. 
# Or, in the cas of an Area threshold of 0.5 km2, the glacier volume covered is 
# 59% and the glacier number reduction is 92% 
rgi <- rgi |> mutate(Vol_km3 = glacierVolume_RGIF(Area))
rgi_cal <- rgi |> 
  dplyr::filter(Area >= 0.1) 

temperature_cal <- hist_obs |> 
  mutate(date = as_date(date)) |> 
  dplyr::select(c(date, starts_with(rgi_cal$RGIId)))

observed_cal <- observed_melt |> 
  rename(RGIID = RGIId, 
         totAbl = Qgl_m3a, 
         year = start) |> 
  mutate(year = year(year)) |> 
  dplyr::filter(RGIID %in% rgi_cal$RGIId)

glaciers_cal <- rgi_cal$RGIId

DDFfit <- NULL
for (id in c(1:length(glaciers_cal))) {
  cat("Fitting temperature index model for glacier", glaciers_cal[id], "id =",
      id, "\n")
  
  tempGA <- ga(type = "real-valued", 
               fitness = glacierDischargeRMSE, 
               temperature = temperature_cal,
               observed = observed_cal,
               index = id,
               lower = c(0.1, 0), upper = c(50, 0), 
               popSize = 500, maxiter = 100, 
               parallel = (parallel::detectCores()-6))
  
  cat("    Result: RMSE = ", -tempGA@fitnessValue, "MF =", tempGA@solution[1], 
      "Tth =", tempGA@solution[2], "\n")
  DDFfit <- rbind(DDFfit, tibble(MF = tempGA@solution[1],
                                Tth = tempGA@solution[2],
                                rsme = tempGA@fitnessValue) |>
                    mutate(ID = glaciers_cal[id]))
save(DDFfit, 
     file = paste0(data_path, "GLACIERS/TIparameters.RData")
}

How to determine the MF parameter for the non-calibrated glaciers? Using the average MF from the calibrated glaciers is one option. The best results (for the present case study) are obtained by modulating the MF according to glacier slope.

Now we re-calculate the glacier melt using the calibrated data set.

# Option to calculate only for a selection of glaciers (modify indices)
indices <- c(1:length(glaciers))  
DDFfit |> dplyr::filter(RGIId %in% glaciers[indices])
## # A tibble: 179 x 9
##       MF   Tth   rsme RGIId          ID               Slope .fitted .resid layer
##    <dbl> <dbl>  <dbl> <chr>          <chr>            <dbl>   <dbl>  <dbl> <dbl>
##  1 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_1  23.7    4.84  -3.93     1
##  2 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_2  23.7    4.84  -3.93     2
##  3 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_3  23.7    4.84  -3.93     3
##  4 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_4  23.7    4.84  -3.93     4
##  5 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_5  23.7    4.84  -3.93     5
##  6 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_6  23.7    4.84  -3.93     6
##  7 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_7  23.7    4.84  -3.93     7
##  8 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_8  23.7    4.84  -3.93     8
##  9 0.916     0 -115.  RGI60-13.06145 RGI60-13.06145_9  23.7    4.84  -3.93     9
## 10 0.547     0  -72.6 RGI60-13.06157 RGI60-13.06157_1  22.8    4.77  -4.23     1
## # ... with 169 more rows
melt_cal <- glacierMelt_TI(temperature = hist_obs |> 
                             dplyr::select(-date),
                           MF = MFvec,
                           threshold_temperature = Tthvec) |> 
  mutate(date = hist_obs$date) |> 
  relocate(date, .before =  where(is.numeric))

melt_a_eb <- melt_cal |> 
  pivot_longer(-date, names_to = "ID", values_to = "Melt") |> 
  mutate(Hyear = hyear(date)) |> 
  group_by(Hyear, ID) |> 
  summarise(date = first(date), 
            Melt = sum(Melt), 
            .lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0, 
                              Melt*(1-error_stats$sMelt)), 
            .ub_Melt = Melt*(1+error_stats$sMelt))|> 
  ungroup() 

melt_a <- melt_a_eb |> 
  separate(ID, into = c("RGIId", "elB"), sep = "_") |> 
  group_by(Hyear, RGIId) |> 
  summarise(date = as_date(first(date)), 
            Melt = mean(Melt), 
            .lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0, 
                              Melt*(1-error_stats$sMelt)), 
            .ub_Melt = Melt*(1+error_stats$sMelt)) |> 
  ungroup() |> 
  dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)) |> 
  left_join(rgi |> dplyr::select(RGIId, Area) |> st_drop_geometry(), 
            by = "RGIId") |> 
  # Assume constant glacier area over the comparison period. 
  mutate(Qgl_m3a = Melt*10^(-3) * Area*10^6, 
         .lb_Qgl_m3a = Qgl_m3a * (1-sqrt(error_stats$sMelt^2+error_stats$sA^2)),   
         .lb_Qgl_m3a = ifelse(.lb_Qgl_m3a < 0, 0, .lb_Qgl_m3a), 
         .ub_Qgl_m3a = Qgl_m3a * (1+sqrt(error_stats$sMelt^2+error_stats$sA^2)))

melt_a_basin <- melt_a |> 
  group_by(Hyear) |> 
  summarise(Qgl_m3a = sum(Qgl_m3a, na.rm = TRUE), 
            .lb_Qgl_m3a = sum(.lb_Qgl_m3a, na.rm = TRUE), 
            .ub_Qgl_m3a = sum(.ub_Qgl_m3a, na.rm = TRUE), 
            date = first(date), 
            Area = sum(Area, na.rm = TRUE), 
            Melt = mean(Melt, na.rm = TRUE)) |> 
  ungroup()

melt_obs_a_basin <- melt_obs_a_all |>  
  mutate(Hyear = hyear(start)) |> 
  group_by(Hyear) |> 
  summarise(Qgl_m3a = sum(Qgl_m3a, na.rm = TRUE), 
            .lb_Qgl_m3a = sum(.lb_Qgl_m3a, na.rm = TRUE), 
            .ub_Qgl_m3a = sum(.ub_Qgl_m3a, na.rm = TRUE), 
            date = first(start), 
            Area = sum(Area_m2, na.rm = TRUE)) |> 
  ungroup()

ggplot() +
  geom_ribbon(data = melt_a_basin,
              aes(date, Qgl_m3a, ymin = .lb_Qgl_m3a, ymax = .ub_Qgl_m3a),
              colour = NA, alpha = 0.1) +
  geom_ribbon(data = melt_obs_a_basin,
              aes(date, Qgl_m3a,
                  ymin = .lb_Qgl_m3a,
                  ymax = .ub_Qgl_m3a, linetype = "obs"),
              size = 0.2, alpha = 0.1) +
  geom_line(data = melt_a_basin,
            aes(date, Qgl_m3a, linetype = "sim")) +
  geom_line(data = melt_obs_a_basin,
            aes(date, Qgl_m3a, linetype = "obs")) +
  labs(x = "Date", y = "Melt [m weq/a]") +
  scale_linetype_manual(name = "Source",
                        values = c("sim" = 1, "obs" = 2)) +
  scale_colour_viridis_d() +
  scale_fill_viridis_d() +
  theme_bw()
Glacier discharge in the Atbashy basin.

Glacier discharge in the Atbashy basin.

Aggregated to the basin level, the simulated glacier discharge fits the observed glacier discharge quite nicely.

comparison <- melt_obs_a_all |> 
  dplyr::filter(RGIId %in% glaciers[indices]) |> 
  mutate(Hyear = hyear(start)) |> 
  rename(Obs = Qgl_m3a, 
         .lobs = .lb_Qgl_m3a, 
         .uobs = .ub_Qgl_m3a) |> 
  left_join(melt_a |> 
              dplyr::filter(RGIId %in% glaciers[indices]) |> 
              rename(Sim = Qgl_m3a, 
                     .lsim = .lb_Qgl_m3a, 
                     .usim = .ub_Qgl_m3a), 
            by = c("Hyear", "RGIId"))

axmin <- min(cbind(comparison$Obs, comparison$Sim), na.rm = TRUE)
axmax <- max(cbind(comparison$Obs, comparison$Sim), na.rm = TRUE)

ggplot(comparison |> drop_na()) +
  geom_abline(slope = 1, intercept = 0) + 
  geom_point(aes(Obs, Sim, colour = RGIId), size = 0.4, show.legend = FALSE) + 
  scale_colour_viridis_d() + 
  xlim(axmin, axmax) + 
  ylim(axmin, axmax) + 
  coord_fixed() + 
  theme_bw()
Observed versus simulated glacier discharge. 11 years of annual data for 131 glaciers in the basin.

Observed versus simulated glacier discharge. 11 years of annual data for 131 glaciers in the basin.

However, when looking at calibration outcomes of individual glaciers, the picture does not look quite as well. For some glaciers, the optimization does not converge to produce good estimates of glacier discharge.

Glacier mass balance

The following figure shows the components of the glacier mass balance for a few glaciers in the Atbashy basin.

glacier_balance <- glacierBalance(melt_a_eb = melt_a_eb, 
                                  rgi_elbands = rgi_elbands)

glacier_balance <- glacier_balance |> 
  mutate(.lb = ifelse(Variable == "A_km2", 
                      ifelse(Value*(1-error_stats$sA)>0, 
                             Value*(1-error_stats$sA), 0), 
                      ifelse(Variable == "V_km3", 
                             ifelse(Value*(1-error_stats$sV)>0, 
                                    Value*(1-error_stats$sV), 0), 
                             ifelse(Variable == "Q_m3a", 
                                    ifelse(Value*(1-error_stats$sQglmelt)>0,
                                           Value*(1-error_stats$sQglmelt), 0), 
                                    ifelse(Variable == "Qimb_m3a", 
                                           Value*(1-error_stats$sImbal), NA)))), 
         .ub = ifelse(Variable == "A_km2", 
                      Value*(1+error_stats$sA), 
                      ifelse(Variable == "V_km3", 
                             Value*(1+error_stats$sV), 
                             ifelse(Variable == "Q_m3a", 
                                    Value*(1+error_stats$sQglmelt), 
                                    ifelse(Variable == "Qimb_m3a", 
                                           Value*(1+error_stats$sImbal), NA)))))

ggplot(glacier_balance |> 
         dplyr::filter(RGIId %in% glaciers[6:7], 
                       Hyear > min(Hyear) & Hyear < max(Hyear), 
                       Variable %in% c("A_km2", "V_km3", "Q_m3a", "Qimb_m3a"))) + 
  geom_ribbon(aes(Hyear, ymin = .lb, ymax = .ub, fill = RGIId), 
              alpha = 0.2, colour = NA) + 
  geom_line(aes(Hyear, Value, colour = RGIId)) + 
  facet_wrap("Variable", scales = "free_y") + 
  scale_colour_viridis_d() + 
  scale_fill_viridis_d() + 
  theme_bw()
Glacier development and discharge.

Glacier development and discharge.

We now have glacier discharge (Q_m3s) and the unsustainable contribution to glacier discharge, the imbalance ablation (Qimb_m3a) which is negative for glacier loss and positive for growing glaciers.

Future glacier melt

scenarios <- names(fut_sim)
fut_melt_d <- NULL
fut_melt_a_eb <- NULL
fut_glacier_balance <- NULL
for (scenario in scenarios) {
  temp_melt_d <- glacierMelt_TI(temperature = fut_sim[[scenario]] |>
                     dplyr::select(-Date),
                   MF = MFvec,
                   threshold_temperature = Tthvec) |> 
      mutate(Scenario = scenario, 
             date = fut_sim[[scenario]]$Date) |>
      relocate(date, .before =  where(is.numeric))  |>
      pivot_longer(-c(date, Scenario), names_to = "ID", 
                   values_to = "Melt") 
  temp_melt <- temp_melt_d |>
      mutate(Hyear = hyear(date)) |> 
      # Sum melt per year
      group_by(Hyear, Scenario, ID) |> 
      summarise(date = first(date), 
                Melt = sum(Melt)) |> 
      ungroup()
  
  fut_melt_d <- rbind(fut_melt_d, temp_melt_d)
  fut_melt_a_eb <- rbind(fut_melt_a_eb, temp_melt)
  
  fut_glacier_balance <- rbind(
    fut_glacier_balance, 
    glacierBalance(melt_a_eb = temp_melt,
                   rgi_elbands = rgi_elbands) |> 
      mutate(Scenario = scenario))

} 
fut_melt_amean <- fut_melt_a_eb |>
  # Mean melt rate per glacier
  # mutate(RGIId = gsub("\\_\\{d}+$", "", ID)) |> 
  # group_by(Hyear, Scenario, RGIId) |> 
  # summarise(Melt = mean(Melt)) |> 
  # Mean melt rate in the catchment
  group_by(Hyear, Scenario) |> 
  summarise(date = as_date(first(date)), 
            Melt_mma = mean(Melt)) |> 
  ungroup() |> 
  separate(Scenario, into = c("GCM", "SSP"), sep = "_") |> 
  mutate(SSP = factor(SSP, levels = c("ssp585", "ssp370", "ssp245", "ssp126"))) |> 
  # dplyr::filter(!str_detect(GCM, "UK")) |>  
  dplyr::filter(!str_detect(GCM, "IPSL"), 
                Hyear < max(Hyear)) |> 
  group_by(Hyear, SSP) |> 
  mutate(lb = min(Melt_mma), 
         ub = max(Melt_mma), 
         mean = mean(Melt_mma)) |> 
  ungroup()

melt_amean <- melt_a |> 
  dplyr::select(Hyear, date, Melt) |> 
  group_by(Hyear) |>  
  dplyr::summarise(date = first(date), 
                   Melt_mma = mean(Melt)) |> 
  ungroup()

ggplot(fut_melt_amean) + 
  geom_ribbon(aes(date, ymin = lb, ymax = ub, fill = SSP), alpha = 0.2) + 
  geom_line(aes(date, mean, colour = SSP)) + 
  geom_line(data = melt_amean, aes(date, Melt_mma)) + 
  # geom_smooth(aes(date, mean, colour = SSP), method = "loess", se = FALSE) + 
  # geom_smooth(data = melt_amean, aes(date, Melt_mma), method = "loess", 
  #             se = FALSE) + 
  scale_colour_viridis_d() + 
  scale_fill_viridis_d() + 
  ylab("Melt [mm/a]") + 
  xlab("Year") + 
  theme_bw()
Future specific glacier melt rates expected in the Atbashy basin.

Future specific glacier melt rates expected in the Atbashy basin.

We see an increase of the expected glacier melt rates. In the two lesser emission scenarios ssp126 and ssp245, glacier melt rates stabilise around mit century and in the last quarter of the century repsectively on a higher level than in recent years while under the higher emission scenarios ssp370 and ssp585 glacier melt continues to rise.

Glacier balance over climate scenarios

Under development.

Uncertainties

The shaded area in the above figure indicates the spread of the GCM ensemble used. The width of the uncertainty band of the model method is at least twice the width of the above uncertainty bands, such that it is near impossible to differentiate between the impact of the different socio-economic scenarios. Dramatic glacier volume reduction cannot be avoided even with the least emissions scenario ssp126.

Peak water

Under development

Glacier volume loss

Under development.

Next steps

Now we have glacier discharge and glacier imbalance ablation. How do we include these in our modelling workflow? Glacier vignette 4 explains how glaciers are treated in the hydraulic/hydrologic model RS Minerve and details a workflow to implement glacier discharge in a model in RS Minerve.

Open Global Glacier Model (OGGM)

The Open Global Glacier Model is a community project which provides open-source glacier modelling tools in python and educational resources. The modeling framework includes comprehensive tutorials and data-pre and post-processing tools. It can be used as an alternative to the glacier modelling work flow presented here.

References

Erasov, N. V. 1968. “Method for Determining of Volume of Mountain Glaciers.” MGI, no. 14: 307–8.
Farinotti, Daniel, Matthias Huss, Johannes J. Fürst, Johannes Landmann, Horst Machguth, Fabien Maussion, and Ankur Pandit. 2019. “A Consensus Estimate for the Ice Thickness Distribution of All Glaciers on Earth.” Nature Geoscience 12 (3): 168–73. https://doi.org/10.1038/s41561-019-0300-3.
Gudmundsson, L. 2016. “Qmap: Statistical Transformations for Post-Processing Climate Model Output.” https://CRAN.R-project.org/package=qmap.
Gudmundsson, L., J. B. Bremnes, J. E. Haugen, and T. Engen-Skaugen. 2012. “Technical Note: Downscaling RCM Precipitation to the Station Scale Using Statistical Transformations – a Comparison of Methods.” Hydrology and Earth System Sciences 16 (9): 3383–90. https://doi.org/10.5194/hess-16-3383-2012.
Hock, Regine. 2003. “Temperature Index Melt Modelling in Mountain Areas.” Journal of Hydrology 282 (1): 104–15. https://doi.org/10.1016/S0022-1694(03)00257-9.
Hugonnet, Romain, Robert McNabb, Etienne Berthier, Brian Menounos, Christopher Nuth, Luc Girod, Daniel Farinotti, et al. 2021. “Accelerated Global Glacier Mass Loss in the Early Twenty-First Century.” Nature 592 (7856): 726–31. https://doi.org/10.1038/s41586-021-03436-z.
Karger, Dirk Nikolaus, Olaf Conrad, Jürgen Böhner, Tobias Kawohl, Holger Kreft, Rodrigo Wilber Soria-Auza, Niklaus E. Zimmermann, H. Peter Linder, and Michael Kessler. 2017. “Climatologies at High Resolution for the Earth’s Land Surface Areas.” Scientific Data 4 (1): 170122. https://doi.org/10.1038/sdata.2017.122.
Krasting, John P., Jasmin G John, Chris Blanton, Colleen McHugh, Serguei Nikonov, Aparna Radhakrishnan, Kristopher Rand, et al. 2018. “NOAA-GFDL GFDL-Esm4 Model Output Prepared for Cmip6 CMIP.” Earth System Grid Federation. https://doi.org/10.22033/ESGF/CMIP6.1407.
Miles, Evan, Michael McCarthy, Amaury Dehecq, Marin Kneib, Stefan Fugger, and Francesca Pellicciotti. 2021. “Health and Sustainability of Glaciers in High Mountain Asia.” Nature Communications 12 (2868): 10. https://doi.org/https://doi.org/10.1038/s41467-021-23073-4.
Scrucca, Luca. 2013. GA : A Package for Genetic Algorithms in r.” Journal of Statistical Software 53 (4). https://doi.org/10.18637/jss.v053.i04.
———. 2017. “On Some Extensions to GA Package: Hybrid Optimisation, Parallelisation and Islands Evolution 9: 20.
Tang, Yongming, Steve Rumbold, Rich Ellis, Douglas Kelley, Jane Mulcahy, Alistair Sellar, Jeremy Walton, and Colin Jones. 2019. “MOHC Ukesm1.0-LL Model Output Prepared for Cmip6 CMIP Esm-piControl.” Earth System Grid Federation. https://doi.org/10.22033/ESGF/CMIP6.5953.
Yukimoto, Seiji, Tsuyoshi Koshiro, Hideaki Kawai, Naga Oshima, Kohei Yoshida, Shogo Urakawa, Hiroyuki Tsujino, et al. 2019. “MRI MRI-Esm2.0 Model Output Prepared for Cmip6 CMIP Historical.” Earth System Grid Federation. https://doi.org/10.22033/ESGF/CMIP6.6842.